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Highly stratified shear layers are rendered unstable even at high stratifications by Holmboe in- 
stabilities when the density stratification is concentrated in a small region of the shear layer. These 
instabilities may cause mixing in highly stratified environments. However these instabilities occur 
in limited bands in the parameter space. We perform Generalized Stability analysis of the two 
dimensional perturbation dynamics of an inviscid Boussinesq stratified shear layer and show that 
Holmboe instabilities at high Richardson numbers can be excited by their adjoints at amplitudes 
that are orders of magnitude larger than by introducing initially the unstable mode itself. We 
also determine the optimal growth that is obtained for parameters for which there is no instability. 
We find that there is potential for large transient growth regardless of whether the background 
flow is exponentially stable or not and that the characteristic structure of the Holmboe instability 
asymptotically emerges as a persistent quasi-mode for parameter values for which the flow is stable. 



I. INTRODUCTION 

The mixing of shear layers and the development of tur- 
bulence is severely impeded when the layer is located in 
regions of large stable stratification. The stratification 
is usually quantified with the non-dimensional Richard- 
son number defined as the ratio of the local Brunt- 
Vaisala frequency iV 2 to the square of the local shear U', 
Ri = N 2 /U' 2 . Large stratification corresponds to large 
Richardson numbers. The significance of the Richardson 
number for the stability of stratified flows has been un- 
derscored with a theorem due to Miles and Howard [1,2] 
which proves that if everywhere in the flow Ri > 1/4, 
the flow is necessarily stable to exponential instability. 
The essential instability that leads primarily to mixing 
in both stratified and unstratified shear layers is the 
Kelvin- Helmholtz (KH) instability [3]. The KH modes 
are eventually stabilized when the local Richardson num- 
ber exceeds 1/4, but if the stratification is concentrated 
in narrow regions within the shear layer the Richardson 
number may locally become smaller than 1 /4 and an in- 
stability can result although the overall stratification is 
very large. Under such conditions, a new branch of in- 
stability of stratified shear layers emerges as shown by 
Holmboe [4], consisting of a pair of propagating waves 
with respect to the center flow, one prograde and one 
retrograde. This new instability branch has been named 
the Holmboe (H) instability and it is physically very in- 
teresting because it persists at all stratifications and can 
lead to mixing in highly stratified shear layers. 

Holmboe instabilities have been reproduced in labora- 
tory experiments [5-11] and have been numerically simu- 
lated [12-19]. At finite amplitude these Holmboe modes 
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equilibrate into propagating waves and they can induce 
mixing in highly stratified environments [12, 17, 20]. 
Consequently, Holmboe instabilities present the intrigu- 
ing possibility that they may be responsible for transi- 
tion to turbulence and mixing in highly stratified flows. 
In that vein, Alexakis [18, 21, 22] has proposed that the 
surmised necessary mixing in order for a thermonuclear 
runway to occur and highly stratified white dwarfs form 
supernova explosions could be effected with Holmboe in- 
stabilities. 

Holmboe [4] also introduced a new method for analyz- 
ing the perplexing instabilities that may arise with the 
introduction of stratification. He went beyond classical 
modal analysis and proposed that a fruitful program for 
predicting and obtaining a clear physical idea of the in- 
stability of stratified shear flows is to simplify the back- 
ground flow to segments of piecewise constant vorticity, 
as in Rayleigh's seminal work [23], and then consider the 
dynamics of the edge waves that are supported at each 
vorticity and density discontinuity. Holmboe showed that 
the flow becomes unstable when the edge waves propa- 
gate with the same phase speed. This method of analysis 
has been extended and used to clarify the detailed mech- 
anism of instability in stratified flows [24-28] and has 
been shown recently by Carpenter et al. [29] to be capa- 
ble of readily assessing the character of instability of non- 
idealized observed flows. Also this edge wave description 
has elucidated the dynamics of other instabilities that 
occur in geophysics [30-33] and astrophysics [34-36]. 

However, the efficacy of Holmboe modes to mix highly 
stratified layers can be questioned because the instability 
occurs only in limited regions of parameter space and the 
modal growth rate of the instability becomes exponen- 
tially small with stratification. Specific estimates of the 
asymptotic behavior of the growth rates and of the un- 
stable band of wavenumbers has been recently obtained 
by Alexakis [22]. Also the introduction of viscosity can 
further reduce the growth rates of the Holmboe insta- 
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bilitics. It is the purpose of this work to go beyond the 
modal analysis and investigate the non-modal stability of 
stratified shear layers that arc Holmboc unstable in order 
to assess the true potential of growth of stratified shear 
layers. We will achieve this using the standard methods 
of generalized stability theory [37, 38]. 

Because shear layers are powerful transient amplifiers 
of perturbation energy, we expect that the Holmboe in- 
stabilities can be excited at enhanced amplitude by com- 
posite non-modal perturbations. Even in the case of an 
infinite constant shear flow which has a continuous spec- 
trum with no inviscid analytic modes and all perturba- 
tions eventually decay algebraically with time, the non- 
modal solutions constructed by Kelvin [39] demonstrate 
that the asymptotic limit is non-uniform and perturba- 
tion energy can transiently exceed any chosen bound in 
the inviscid limit; the same is true for bounded Couette 
flow as shown by Orr [40]. The Kelvin-Orr solutions can 
be extended to stratified flows [41, 42] and these solutions 
can produce transient perturbation growth that can also 
exceed in the inviscid limit any bound [43]. Specifically 
it can be shown that in constant shear flow and large 
Richardson numbers the perturbation energy amplifica- 
tion is the square root of that achieved in the absence 
of stratification. The same robust growth resulting from 
the continuous spectrum is expected to persist in all shear 
layers as the dynamics of the Orr waves are universal and 
do not depend on the details of the background shear 
flow [44]. 

Furthermore in shear layers of finite size the vortic- 
ity and density edge waves that are concentrated in re- 
gions of vorticity and density discontinuities can interact 
and produce transient perturbation growth [45] . Also the 
transiently growing Kelvin-Orr vorticity structures may 
deposit their energy to the edge waves [33] or to prop- 
agating gravity waves [46-48]. In general the continu- 
ous spectrum can excite at high amplitude the unstable 
or even stable analytic modes [49, 50] and maintain un- 
der continuous excitation high levels of variance in stable 
flows [51, 52]. 

In this work we will consider the generalized stability 
of a Boussinesq stratified shear layer. The Boussinesq 
approximation makes inaccurate predictions of the dis- 
persion relation of the perturbations in the shear layers 
when the stratification is locally very large [53]. How- 
ever we have found that the resulting optimal growths 
are not sensitive to this approximation and in this pa- 
per in order to make contact with previous work we will 
only present results obtained using the Boussinesq ap- 
proximation. Also following Hazel [54] and the work of 
Smyth and Peltier [55, 56] we will consider a continuous 
version of Holmboc's velocity profile with the character- 
istic density stratification embedded in the shear layer. 
Hazel and Smyth and Peltier had analyzed the charac- 
ter of the bifurcations that occur in the instability of 
the continuous profile as the Richardson number at the 
center of the shear layer increases. Alexakis [18, 21, 22] 
extended the analysis to even higher stratifications and 



was able to predict theoretically and show numerically 
that there are higher branches of Holmboe instabilities 
in these continuous profiles. Because of the geophysical 
and astrophysical interest we will concentrate our anal- 
ysis of non-modal perturbation growth of shear layers at 
very high Richardson numbers which possess sparse is- 
lands of exponential instability of low growth rate. 



II. 



FORMULATION 



We will study the stability of an inviscid, unidirec- 
tional, stratified, step like channel flow to two dimen- 
sional incompressible perturbations. The mean flow is 
characterized by a velocity jump AU over a vertical scale 
ho and the associated statically stable mean density is 
characterized by a density jump Ap. The linearized 
nondimensional perturbation equations about the mean 
flow Uq(z) and mean density po(z) in the Boussinesq ap- 
proximation are: 

(d t + U d x ) u = -Uq w - d x p (la) 

(d t + U Q d x )w = -Jg-d zP (lb) 

(d t + U d x ) g = -p' w (lc) 

d x u + d z w = (Id) 

The density of the fluid has been decomposed as: p = 
Pm + Po(z) + g(x, z, t), where p m is the mean background 
density, po{z) is the mean density variation and g(x, z, t) 
is the perturbation density; and furthermore, according 
to the Boussinesq approximation, we assume, \po\ *C p m 
and \g\ <C p m - We denote with u and w the perturbation 
velocities in the strcamwise (x) and vertical (z) direc- 
tion and p is the perturbation reduced pressure. Dif- 
ferentiation with respect to z is denoted with a dash. 
The perturbation equations have been made nondimen- 
sional, choosing ho as the length scale, AU as the velocity 
scale, and Ap as the scale for the density. Time has been 
scaled with ho/AU, and pressure has been scaled with 
p m (AU) . The local Richardson number is 



with 



Ri 



J = 



-J 



Po 
II' 2 



9 Ap hp 
p m AU 2 



(2) 



(3) 



J provides a measure of the bulk Richardson number of 
the shear layer. We will assume that the flow is in a chan- 
nel of length L = 2ho and impose at the channel bound- 
aries zero vertical velocity, w, and zero perturbation den- 
sity, g. The length of the channel has been selected so 
that boundaries do not affect in an important way the 
results that we report in this paper. We have checked 
by doubling the channel size that the eigenspectrum and 
the perturbation transient growth are converged. Similar 
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insensitivity to the location of the channel walls has been 
previously reported by Alexakis [18]. 

The perturbation equations can be written compactly 
by introducing a streamfunction tp so that the velocities 
are u = d z ip and w = —d x ip, and the vertical displace- 
ment, defined through the relation: 



(d t + U d x ) r\ = w , 
can replace the perturbation density as: 
Q = -Po V ■ 



(4) 



(5) 



The perturbation state will be determined by the fields ip 
and r\. Because of homogeneity in the strcamwise direc- 
tion, the perturbation state can be decomposed in Fourier 
components in x: 



x(z,t)e ik * = $(z,t),fj(z,t) 



T e ikx 



(6) 



and the perturbation dynamics for the Fourier ampli- 
tudes for each wavenumber k can be written compactly 
in the form: 



dx 
~&t 



with the operator A given by: 
A = -ik 



A- 1 (U A - US) 
1 



JA 



u 



(7) 



(8) 



With A = &l z — k 2 we denote the two dimensional Lapla- 
cian for the k streamwise wavenumber. A -1 is the in- 
verse Laplacian, the inversion being rendered possible 
and unique by imposing the boundary conditions on the 
upper and lower boundaries. The perturbation equation 
(7) is equivalent to the time dependent Taylor- Goldstein 
equation: 



(d t + ikU ) 2 - ikUo (d t + ikU ) V> + fc 2 Jp'$ = 



(9) 



III. KELVIN-HELMHOLTZ AND HOLMBOE 
INSTABILITIES OF A SHEAR LAYER 

Following Hazel [54], we consider the following back- 
ground mean velocity and density: 



Uo(z) = ^tanh(2z) , p (z) 



- tanh (2Rz) . (10) 



The mean velocity and density profiles are shown in 
FIG. 1 for R = 3. The parameter R is the ratio of the 
width of the shear layer to the width of the region over 
which the density jump occurs. The local Richardson 
number is given by, 



RAO) = RJ 



scch 2 (2Rz) 
sech 4 (2z) 



(11) 



The Richardson number at the center is related to the 
bulk Richardson number through the relation: 



Ri(0) = RJ 



(12) 



When R = 1 the local Richardson number assumes its 
minimum value at the center, Ri(0) = RJ, and the shear 
layer is unstable only if Ri(0) < 1/4. The instability that 
results is a Kelvin- Hclmholtz (KH) type instability with a 
short wave cutoff and zero phase velocity. For R > 2 the 
Richardson number decays exponentially to zero. This 
exponential decay allows the local Richardson number to 
be in regions smaller than 1/4 and the shear layer may 
become unstable although both the bulk stratification J 
and Ri(0) may be very large. The instabilities that occur 
under these conditions are the Holmboe instabilities (H). 
The variation of the Richardson number with height is 
shown in FIG. 2 for R — 3, a case that supports Holmboe 
instabilities and will be treated in detail in this paper, 
and for R = \/2 which can support only KH instability 
when the Richardson number at the center is smaller than 
1/4. 

In order to determine the modal stability of the veloc- 
ity and density profiles (10) we assume modal solutions 
of (7) of the form x(z,t) = [ip c (z),fi c (z)}' T e~ lkct , with 
— ike an eigenvalue of the evolution operator A in (8). 
The modal growth rate of the perturbation is kci where 
Cj = S(c) and the flow is exponentially unstable to per- 
turbations of wavenumber k if a > 0. The real part, c r , 
of the eigenvalue c(k), gives the phase speed of the per- 
turbation. Detailed analysis of the eigenspectrum of (10) 
and the bifurcation from KH instabilities to H instabili- 
ties has been already performed [18, 21, 22, 54-56]. Here 
we review the basic stability characteristics for the case 
of R = 3 that admits both KH and H instability modes. 

The transition from KH instability to H instability is 
shown in FIG. 3 where we plot the growth rate and phase 
speed as a function of the Richardson number at the cen- 
ter of the shear zone Ri(0) for k = 0.3. For Ri(0) < 0.242, 
there is a single unstable Kclvin-Hclmholtz instability 
with phase speed c r = 0. At Ri(0) = 0.242 a second 
KH mode becomes unstable, also with c r = 0. The KH 
modes coalesce at Ri(0) = 0.258 to form a propagating 
pair of unstable Holmboe modes with equal and opposite 
phase speeds and the same growth rate. In this case we 
see two branches of KH instabilities existing even when 
the center Richardson number exceeds 1/4 (the Richard- 
son number assumes values smaller than 1/4 way from 
the center because the parameter is R — 3, see FIG. 2); 
a similar case can be found in Smyth and Peltier [55]. 

Numerical results were obtained by discretizing the 
channel and all differential operators using second order 
centered differences and incorporating the appropriate 
boundary conditions. Equation (7) then becomes a ma- 
trix equation in which the state becomes a column vec- 
tor. We used a staggered grid, that is, we evaluated ip c (z) 
at N interior points in the vertical and fj c (z) at N + 1 
points located halfway between the collocation points of 
the streamfunction. The eigenspectrum is obtained by 
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cigcnanalysis of operator A. Calculations, unless other- 
wise specified, were performed with TV = 1001 in the do- 
main — 2 < z < 2. With this discretization the finite time 
evolution of the dynamics is well resolved when there is 
no instability up to a time of t = 0(l/(akAz)) where a 
is the typical shear [43]. However, in order to numerically 
resolve the curves of zero modal growth of the Holmboc 
modes of instability wc had to include numerical diffu- 
sion in both the momentum and density equations with 
coefficient of diffusion v = (Az/2) 2 in the momentum 
equation and u/9 in the density equation for the case 
R = 3 (Az is the grid spacing). Very accurate determi- 
nation of the neutral curve in the inviscid limit can be 
obtained, using the shooting method as in Alexakis [22]. 

Growth rate contours as a function of stratification, J 
or Ri(0), and wavenumber, k are shown in FIG. 4 for 
R — 3. The region of KH instability is concentrated for 
small Richardson numbers and it is not discernible in this 
plot. However the bands of instability corresponding to 
the various Holmboe modes of instability can be clearly 
seen. Note that because of the inclusion of diffusion the 
instability bands do not extend to infinity as they do in 
the inviscid limit [4, 21, 24, 57]. 

The first band in FIG. 4 is the first Holmboe mode of 
instability, HI, (for J > 0.0833), the second is the second 
Holmboe instability mode, H2, (for J > 3.5) and we can 
also partly see the growth rate of the third Holmboe in- 
stability mode, H3 (for J > 10). Estimates of the growth 
rates of these higher modes and regions of instability can 
be found in Alexakis [22] who shows that the growth rate 
of the Holmboc instabilities is decreasing exponentially 
with stratification and as the Holmboe modes approach 
neutrality the critical layer of the instabilities, that is the 
height z c for which Uq(z c ) — c r , tends to the height that 
corresponds to the maximum/minimum velocity ±0.5 of 
the shear zone [22] . We obtain a sense of the ordering of 
the growth rates of the KH and H instabilities in FIG. 5. 
In that figure we plot the modal growth rates kci as a 
function of wavenumbcr k for R = 3 and center Richard- 
son numbers Ri(0) = 0.06, 0.65, 12 for which we have 
respectively KH mode instability, HI and H2 instabili- 
ties. 

Typical structure of the KH instability mode is shown 
in FIG. 6 for k = 1, R = 3 and Ri(0) = 0.06. For these 
parameters the coalescence of the KH modes and the 
emergence of the first Holmboe mode of instability occurs 
at Ri(0) = 0.182. The Holmboe modes come in prograde 
and retrograde pairs. The prograde and retrograde un- 
stable Holmboe modes, HI, for Ri(0) = 0.65 are shown 
in FIG. 7. Because they are symmetric with respect to 
z = 0, wc will subsequently plot only the prograde branch 
of the H mode. Higher Holmboe modes have larger ver- 
tical wavenumbers, the case of H2 is shown in FIG. 8. 
In the idealized piecewise constant problem studied by 
Holmboe (see Appendix A) the prograde Holmboe insta- 
bility arises from resonant interaction between the edge 
waves at the density discontinuity at the center with the 
edge wave at the upper vorticity discontinuity, and with 



the edge wave at the lower vorticity discontinuity for the 
retrograde instability [24]. The same type of interaction 
gives rise to the Holmboc instability modes when the dis- 
continuities of the background are smoothed as for the 
background given by Eq. (10). Because of the delta func- 
tion nature of the vorticity and stratification gradient the 
edge waves in Holmboe's idealized profile do not have in- 
ternal structure and as a result there is a single pair of 
unstable Holmboc instabilities that can arise from the 
interaction. The smoothed profile allows the edge wave 
to obtain vertical structure within the regions of large 
vorticity and stratification gradient (cf. respectively to 
the left and right panels of FIG. 7 and FIG. 8) and the 
H2, H3, etc Holmboe instabilities emerge, with the higher 
modes associated with higher vertical wavenumber struc- 
ture. 



IV. OPTIMAL EXCITATION OF 
KELVIN-HELMHOLTZ AND HOLMBOE MODES 

The perturbation dynamics in a stratified shear flow 
are non-normal and, as we will show in this section, the 
unstable modes emerge through excitation of their ad- 
joint. In order to determine the growth that results by 
introducing initially the adjoint of the unstable mode we 
derive the adjoint perturbation dynamics in the energy 
norm. For a perturbation of wavenumber k the total 
average kinetic energy over a wavelength in the region 
[zi,z 2 ] is: 



22 



and the potential energy is: 



= | idzJ(-p' Q )fffi, 



(13) 



(14) 



The total energy of a perturbation x = [i/j, r)] T is thus 
given by the inner product: 



1 fZ2 

■ ' dz 



1*11 = (*,*) = 



V> (-AVi) + J(-p' )fj*v 



(15) 

The adjoint operator in this inner product is defined 
as the operator that satisfies: 



f Q ,g = f Q ,Ag 



(16) 



for any two states f Q = [tpccVa)' 1 and g = [tjj,fj\ T . For a 
discussion of the adjoint equation in fluid mechanics see 
Farrell [58], Hill [59], Farrell and Ioannou [37], Chomaz 
ct. al. [60] and Schmid and Hcnningson [38]. 

From (16) the adjoint operator of (8) in the energy 
inner product is easily shown to be: 



A" 1 (C/ A + 2l^D) 



JA- 1 ^ 
U 



(17) 
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with D = d/dz denoting differentiation with respect to 
z, and the adjoint evolution equation is 



dig 

at 



(18) 



with the adjoint boundary conditions x a (zj,£) = at the 
channel walls z%. This equation produces the adjoint of 
the time-dependent Taylor- Goldstein equation (9): 

(d t - ikUof A$ a -2ikUo (d t - ikU ) Bijj a — k 2 Jp' 'ijj a = 

(19) 

The eigenvalues of the adjoint operator A^ are the com- 
plex conjugate of the eigenvalues of A and in this way 
we can establish a correspondence between the modes of 
the adjoint A^ and the modes of A. Further, the mode of 
the adjoint with eigenvalue ike*, x QC , and the mode of 
A with eigenvalue —ike', x c <, satisfy the biorthogonality 
relation: 



(x QC ,x c ^ =0 for c^c' 



(20) 



From (19) it can be verified that for any analytic mode of 
the Taylor- Goldstein equation (9) with eigenvalue — ike 
(with ci > or with c r outside the range of the back- 
ground flow) the mode of A^ with eigenvalue ike* is: 



x 



4'* c 



Uo-c* 
V* 



(21) 



From the biorthogonality relation (20) arises the phys- 
ical importance of the modes of the adjoint: of all pertur- 
bations of unit energy the largest projection on a given 
mode is achieved by the adjoint of the mode. Indeed, a 
state written as a superposition of the modes of A: 



will have coefficients a c given by: 

(Xq c j x) 



(Xq c j X c ) 



(22) 



(23) 



Because | (x ac ,x) | < 1, if the mode and the state x are 
normalized, with equality when x and x Q c arc multiple 
to each other, the coefficients satisfy the inequality: 



1 



ip^-a Ci ^c) I 



(24) 



and the maximum projection on a given mode is achieved 
by the adjoint mode, i.e. by choosing x = x QC . The 
adjoint perturbation excites the mode at an energy which 
is a factor 



1 



{^-a c j X c j | z 



(25) 



greater than an initial condition consisting of the mode 
itself. For highly non- normal systems |(x QC ,x c )| <C land 
as a result this amplification may be very large, implying 
that the emergence of the mode in these systems is mainly 
due to the excitation of the adjoint. Although this is not 
a new result, it has not been noted in previous studies 
concerning Holmboe instabilities. 

In FIG. 9 we plot the structure of the adjoint in the 
energy inner product of the unstable KH mode shown in 
FIG. 6 with k = 1, Ri(0) = 0.06 and R = 3. An initial 
condition consisting of the adjoint of the most unstable 
mode excites the mode with energy 7.5 times greater than 
an initial condition consisting of the unstable mode itself. 
In FIGs. 10 and 11 we plot the adjoints of the two un- 
stable Holmboe HI and H2 modes, shown in FIGs. 7 and 
8, both with the same k = 1, R = 3 and stratifications 
Ri(0) = 0.65 and Ri(0) = 12 respectively. The adjoints 
of the HI and H2 modes are concentrated at the critical 
layer of the modes which is located way from the center 
of the shear layer and located in a region at the wings of 
the shear layer where the local Richardson number is less 
than 1/4. An initial condition consisting of the adjoint 
of the HI mode excites the most unstable mode with en- 
ergy 69 greater than an initial condition consisting of the 
unstable mode itself. This energy amplification reaches 
7921 for the excitation of the H2 mode by its adjoint. 

Initial conditions in the form of the adjoints of the 
most unstable mode evolve changing form and eventu- 
ally assume the corresponding modal form. During their 
evolution, the adjoint perturbation extracts energy from 
the mean flow which is eventually deposited to the mode 
itself exciting it in this way at high amplitude. The time 
development of the energy of a unit energy initial condi- 
tion in the form of the adjoints of the HI and H2 mode 
is plotted in FIG. 12 demonstrating the increased exci- 
tation of the corresponding modes. This demonstrates 
that the unstable Holmboe modes at high Richardson 
numbers arise primarily due to excitation by their ad- 
joint and the modal growth underestimates the growth 
of the instabilities. 



V. OPTIMAL GROWTH OF PERTURBATIONS 
AND THE EMERGENCE OF THE HOLMBOE 
QUASI-MODE 

In the previous section we have demonstrated that for 
large Richardson numbers the adjoint of a mode can ex- 
cite the modes at much higher amplitude. In this section 
we investigate the optimal growth of perturbations. The 
optimal growth[37, 38, 58] of perturbations at time t in 
the energy metric is obtained by calculating the norm of 
the propagator e A,Mt where A M = M 1 / 2 AIM~ 1 / 2 and IM is 
the energy metric defined as 



Az 



-(D 2 




-Jp'o 



(26) 
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so that total perturbation energy is given by E = xt|Hx. 
The optimal growth is given by the largest singular value, 
<T mM , of e A|hlt and determines the largest perturbation 
growth that can be achieved at time t. The optimal 
perturbation, the initial perturbation that produces this 
growth, is x opt = IM _1 / 2 v, where v the right singular 
vector with singular value <7 ma x- 

We calculate the optimal growth for two indicative op- 
timizing times t = 100 and t = 600 and plot the finite 
time Lyapunov exponent In (cr max (t))/t associated with 
the optimal perturbations as a function of wavenumber 
k and stratification J. As t — > oo the exponent tends 
to the modal growth rate ln(cr max (4))/£ — s- fccj. We saw 
that the modal growth rate (cf. FIG. 4) is non zero only 
in narrow bands of parameter space. In contrast, the 
growth rate associated with the optimal perturbations 
covers all parameter space. For example, for optimiz- 
ing time t = 100 the finite Lyapunov exponent, shown 
in FIG. 13, is almost constant for large Richardson num- 
bers and the equivalent growth rates for large Richardson 
numbers are at least an order of magnitude larger than 
the modal growth rates for all values of the parameters. 
Furthermore the growth rates do not reveal the underly- 
ing bands of exponential instability. The optimal growth 
is robust even for optimizing time t = 600. Contours of 
the finite Lyapunov exponents for this optimizing time 
(FIG. 14) show that the optimal growth rates continue to 
be at least an order of magnitude larger than the modal 
growth rate for all parameter values. Moreover the op- 
timal growth rates are almost constant as a function of 
Richardson number for large Richardson numbers. This 
can be seen in FIG. 15 where we compare the optimal 
growth rate for optimizing times t = 50, 200, 600 with 
the modal growth rate as a function of the stratification 
at the center Ri(0) for k = 1. While the modal growth 
rate is only substantial for the HI branch of instability 
the optimal growth rates produce sustainable growth at 
all stratifications and is reduced by a factor of less than 
2 when the center Richardson number increases tenfold 
from 5 to 50. 

The structure of the optimal perturbations for the 
larger optimizing times is very close to the structure of 
the adjoint of the modes, even for wavenumbers for which 
there is no Holmboe instability. This reveals that excita- 
tion of these flows at high Richardson numbers will lead 
to the emergence of propagating perturbations which are 
close in structure to Holmboe waves. To be specific, con- 
sider central stratification Ri(0) = 5 and excitation of 
the shear layer with various wavenumbers k. In FIG. 16 
we plot the growth rate as a function of k. There is 
a band of instability producing HI waves with phase 
speeds close to c r « ±0.5 and when the flow is excited at 
these wavenumbers pairs of propagating waves emerge. 
This can be seen in the Hovmoller diagram FIG. 17(b) 
in which contours of the logarithm of the positive real 
part of 3t(^(z, t)e tkx ) for k = 3.5 and fixed z = 0.7 are 
plotted in the (x, t) plane for an initial condition in the 
form of the adjoint of the unstable HI mode. The char- 



acteristics of the prograde propagating Holmboe wave 
emerge from the start. In FIG. 17(a) we plot the cor- 
responding Hovmoller diagram when the same adjoint 
is introduced as initial condition and evolved this time 
with the dynamics for k = 1.75. Although at k = 1.75 
there is no instability, a propagating structure emerges 
with the characteristics of the unstable Holmboe wave. 
The same propagating structure also emerges at other 
wavenumbers k. Therefore the Holmboe wave is a robust 
dynamic entity that forms a quasi-mode at wavenumbers 
which do not support an unstable Holmboe wave. This 
quasi-mode is the manifestation of the near resonant edge 
wave structures that lead to the Holmboe instability at 
resonance[4, 24-27, 29]. Similar quasi-mode behavior is 
seen in the evolution of the edge waves in the Holmboe 
profile (cf. Appendix A) for wavenumbers k for which 
there is no instability. When k is close to the wavenum- 
ber for which there is instability the edge wave complex 
periodically amplifies while propagating at a phase speed 
close to the phase speed of the unstable Holmboe modes. 
What is surprising here is that these propagating quasi- 
modes can be excited at such high amplitude and that 
their growth is so persistent. This is shown in FIG. 18 
where the adjoint of the Homboc instability for k = 3.5 
excites at high amplitude quasi-modes for the wavenum- 
bers k = 1, 2, 2.5, 4 for which the flow is stable. The 
same highly amplified quasi-mode would emerge if we 
initialize the flow with a large time optimal. 

The structure of the quasi-mode can not be ascribed 
to the structure of any single mode of the operator. The 
quasi-mode is the superposition of a multitude of con- 
tinuum spectrum modes. The quasi-mode can propagate 
as a coherent entity at a single phase speed because the 
distribution of the amplitude of the coefficients of the 
modal expansion is sharply peaked at the phase speed of 
the propagation of the quasi-mode (the amplitude of the 
coefficients is time invariant because all modes have zero 
growth). For example, the distribution of the amplitude 
of coefficients of the modal expansion of the quasi-mode 
for k = 1.75 in FIG. 19, is concentrated at the phase 
speed c r « ±0.4278 which is exactly the phase speed 
that emerges in FIG. 17(a). 

The quasi-mode can be identified from the frequency 
response of the perturbation dynamics to harmonic forc- 
ing. Because the operators are either neutral or unstable 
in order to obtain steady harmonic response we introduce 
a linear damping in the dynamics, that docs not affect 
the eigenstructures, and consider the frequency response 
of the stable operators Am = Am — 0.02 1, where 1 is the 
identity. The steady state harmonic response, ip(z,w), 
where ip(z,t) = ip(z, w)e lw *, to harmonic forcing Fe iU}t is 
then given by 

$(z,w) = R fc (w)F, (27) 
where Rfc(w) is the resolvent given by 

R fc (w) = (^1-Am)- 1 . (28) 
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The maximum possible response to harmonic forcing is 
then given by the square of the 2-norm of the resolvent, 
||IRfe(a;)|| 2 , which is equal to the square of its largest sin- 
gular value. This determines the maximum energy that 
can be achieved at frequency u> by unit energy harmonic 
excitation. The frequency response shown in FIG. 20 
demonstrates the concentration of the response at the 
phase speeds of the emergent quasi-modes which is very 
close to the phase speed of the Holmboe waves at k = 3.5. 
It demonstrates also the importance of the non-normal 
interactions in the excitation at high amplitude of the 
quasi-modes. In the same plot we plot the equivalent 
normal frequency response 



1 



max 



ikcj\ 2 



(29) 



where ikcj is the j-th eigenvalue of Am that would have 
resulted if the eigenvectors were orthogonal. The differ- 
ence between the responses reflects the excess energy that 
is maintained by the system against friction because of 
the non-orthogonality of the eigenmodes[52, 61]. 



VI. CONCLUSIONS 

Highly stratified shear layers are susceptible to Holm- 
boe instabilities which can lead to mixing. We have 
demonstrated in this work that especially at large 
Richardson numbers the adjoint of the weakly unstable 
modes can excite, through potent transfer of energy from 
the mean flow, the unstable modes at high amplitude. 
Further we found, that even in regions of parameter space 
where the flow is neutral, optimal perturbations can grow 
strongly and excite at large amplitude long-lived prop- 
agating quasi-modes. We have demonstrated that the 
modal growth substantially underestimates the growth 
potential of perturbations in such highly stratified shear 
layers. 
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Appendix A: Non-modal growth produced by the 
edge waves in the idealized Holmboe background 
state 



We consider perturbations in Holmboe's[4] idealized 
piecewise linear velocity profile Uq(z) = z for \z\ < 1/2 
and Uo(z) = z/(2\z\) for \z\ > 1/2 and mean density 
po(z) = —z/\z\ in an infinite domain. This profile has 
two vorticity discontinuities at z = ±1/2 and a density 
discontinuity at z = 0; each vorticity discontinuity sup- 
ports an edge wave and the density discontinuity sup- 
ports a pair of edge waves one prograde and the other 
retrograde [4, 24, 27, 28]. Consider perturbations with 
streamfunction ip(z, t)e lkx . A time dependent solution 
to the perturbation equations can be obtained[4, 24] by 
introducing solutions of the form, 

ip 1 {z,t) = A 1 (t)e- k( - z - 1/2 \ for z > 1/2 (Ala) 

$ 2 CM) = A 2 (t)e~ kz +B 2 (t)e kz , for < z < 1/2 

(Alb) 

$ 3 (z,t) = A 3 (t)e~ kz + B 3 (t)e kz , for - 1/2 < z < 

(Ale) 

Mz, t) = B 4 (t)e k( - z+1 W , for < z < 1/2 (Aid) 

which reduce them, by imposing the usual continu- 
ity conditions at the discontinuity interfaces, to a set 
of ordinary differential equations for the evolution of 
the amplitudes A and B. We choose as a state vari- 
able u(t) = [A 1 (t),C 2 (t),B 4 (t), Vo (t)] T , where m {t) = 
f]{z = 0,t) is the displacement at z = and C 2 (t) = 
(A 2 (i) + B 2 (t)) e fc / 2 . The perturbation equations are 
thus reduced to: 
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du 



(A2) 



with A = IK 1 L with the matrices IK and L given by: 



-e 
2 





-2 1 



(l + e fc (fc-l))/2 




2 

-e* 




-k/2 


k/2 
-ke~ k l 2 



(l + e fc (fc-l))/2 



J ( e fc/2 _ e -fc/2) 



(A3a) 



(A3b) 
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These equations determine the time evolution of the per- 
turbation structure as determined by the interaction of 
the four edge waves. It determines fully the modal sta- 
bility properties of the background state and the non- 
normal dynamics that derive from the interaction of 
the edge waves. In this formulation we have neglected 
the continuum spectrum. A contour plot of the result- 
ing growth rate is shown in FIG. 21 as a function of 
wavcnumber, k, and stratification parameters, J. The 



KH branch of instability is at low k and J and the nar- 
rowing band of the Holmboe mode of instability at larger 
k and J. 

In order to study the non-normal dynamics associated 
with the edge waves we must introduce the correspond- 
ing to (15) energy metric. The potential energy, since 
p' = —S(z), is P — j J 1 770 1 2 - The kinetic energy can be 
written as: 



T =- x ( 



z=l/2 



+ - x ( V3 D^3 



z=0 



+ ip 2 



1p3 Vtp 3 



--1/2 



1p 2 D^2 



= -1/2 



+ 1p4 D'04 



= -1/2 



(A4) 



We thus derive that the perturbation energy of the state 
u(t) is given by 



E = u^u 
with the metric, HI, given by 



where HIt is: 



1 

M = 4 



2k 







(e k - 1) 



-1 1 







-1 
+ e 
-1 



(A5) 



(A6) 



(A7) 



The non-modal growth associated by the dynamics of 
the edge waves is obtained by calculating the optimal 



energy growth that can be achieved at time t which is 
given by: 



E opt {t) = ||exp(A M i)|| 5 



(A8) 



where A M = M 1/2 AM ~ 1/2 . It can be shown that the edge 
waves can lead to substantial growth for parameter values 
for which there is no instability but are close to the sta- 
bility boundary. For example consider the large value of 
stratification J — 1.5. The flow supports unstable waves 
only for 4.2528 < k < 5.142 and the maximum growth 
rate is 0.064. The optimal energy growth produced by the 
edge waves for k = 4.2528, 4.24, 4 is shown in FIG. 22 to 
be substantial as the stability boundary is approached. 
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FIG. 2: Richardson number as a function of height for the velocity and density profiles in FIG. 1 for J = 0.2 and R = 3 (solid). 
The Richardson number assumes values smaller than 1/4 away from the center and the flow may be exponentially unstable. 
Also shown is the case of J = 0.2, R=V2 (dashed). In this case the Richardson number is everywhere greater than 1/4 and 
the velocity profile is by necessity exponentially stable. The threshold value Ri = 1/4 (dash-dot) is also indicated. 



11 




FIG. 3: (Color online) Bifurcation properties of the unstable modes as a function of the Richardson number at the center of the 
shear zone Ri(0) = RJ for the mean state (10) for R = 3 and perturbations with k = 0.3. (a) the growth rate fee, (solid); (b) 
the associated phase speed c r (dashed-dot). For Ri(0) < 0.242 there is a single Kelvin-Helmholtz (KH) mode with zero phase 
speed. For 0.242 < Ri(0) < 0.258 there are two KH modes with zero phase speed. The two KH modes coalesce to produce two 
unstable Holmboe (H) modes with equal and opposite phase speeds. 




FIG. 4: (Color online) Contours of modal growth rate kci of the perturbation operator A as a function of wavenumber, k, and 
the bulk Richardson number, J, for R — 3. Because of the inclusion of diffusion the instability bands do not extend to infinity 
as the do in the inviscid limit. The corresponding value of the Richardson number at the origin Ri(0) = RJ is also recorded at 
the ordinate axis on the right. For this value of R there are KH unstable modes with phase speed c r = for small values of J 
(barely discernible in this plot) and three unstable bands of pairs of H instabilities with prograde and retrograde phase speeds 
for higher J. 
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FIG. 5: (Color online) Growth rates kd as a function of wavenumber k for: a KH instability for Ri(0) = 0.06 (dash-dot), an 
HI mode of instability for Ri(0) = 0.65 (solid) and an H2 mode of instability for Ri(0) = 12 (dashed). The growth rates of H2 
mode have been multiplied by 10. 



Ri(0) = 0.06 

(a) c = 0. 0000 + 0. 1482i (b) 

T=0.87 V=0.13 




FIG. 6: (Color online) Structure of the streamfunction (panel (a)) and density (panel (b)) perturbation fields of an unstable 
KH mode with growth rate kci — 0.1482 and phase speed c r = 0. Critical layer at z c — (dashed). For perturbations with 
k = 1, and base flow with R — 3 and Ri(0) = 0.06. 
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Ri(0) = 0.65 

(a) c = 0. 2423 + 0. 0331i (b) 

T=0.49 V=0.51 




Ri(0) = 0.65 

(c) c = -0. 2423 + 0. 0331i (d) 

T=0.49 V=0.51 




FIG. 7: (Color online) Structure of the streamfunction and density perturbation fields of unstable HI modes with growth rate 
ka = 0.0331 and phase speeds c r = ±0.2423. Panels (a) and (b): streamfunction and density contours for the prograde HI 
mode. Panels (c) and (d): for the retrograde HI mode. Critical layer at z c = ±0.264 (dashed). For perturbations with k = 1, 
and base flow with R = 3 and Ri(0) = 0.65. 
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FIG. 8: (Color online) Structure of the streamfunction (panel (a)) and density (panel (b)) perturbation fields of unstable 
prograde H2 mode with growth rate ka = 0.0006 and phase speed c r = 0.4425. Critical layer at z c = 0.698 (dashed). For 
perturbations with k = 1, and base flow with R = 3 and Ri(0) = 12.00. 




FIG. 9: (Color online) Structure of the streamfunction (panel (a)) and density (panel (b)) perturbation fields of the adjoint of 
the unstable KH mode of FIG. 6. Parameters as in FIG. 6. 
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FIG. 10: (Color online) Structure of the streamfunction and density perturbation fields of the adjoint of the unstable prograde 
HI mode of FIG. 7. The adjoint is centered at the steering level, z c = 0.264, of the mode. 



Ri(0) = 12.00 
(a) c = 0. 4425 + 0. 0006i (b) 

T=0.49 x io" 3 V=0.51 




FIG. 11: (Color online) Structure of the streamfunction and density perturbation fields of the adjoint of the unstable prograde 
H2 mode of FIG. 8. The adjoint is centered at the steering level, z c = 0.698, of the mode. 
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FIG. 12: (Color online) Time evolution of the perturbation energy when the adjoint of the HI mode shown in FIG. 10 and the 
H2 mode shown in FIG. 11 are introduced at t — (solid); also shown is the energy evolution if the HI mode and H2 mode 
are introduced at t = (dashed). The adjoint excites the corresponding modes optimally. The adjoint of the HI mode excites 
the mode at an amplitude a factor 69 greater in energy than an initial condition consisting of the HI mode itself. The H2 
adjoint excites the mode at an amplitude a factor 7921 greater in energy. Parameters are k = 1, R = 3 and for the HI mode 
the stratification is Ri(0) = 0.65, while for the H2 mode Ri(0) = 12. 




FIG. 13: (Color online) Optimal growth for t = 100. Contours of finite time Lyapunov exponent In (<r ma x(t))/t for t = 100 as a 
function of wavenumber, k, and the bulk Richardson number, J, for R = 3. The corresponding value of the Richardson number 
at the origin Ri(0) = RJ is indicated at the ordinate axis on the right. 
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FIG. 14: (Color online) Optimal growth for t = 600. Contours of finite time Lyapunov exponent In (o" m ax(t))/t for t = 600 as a 
function of wavenumber, k, and the bulk Richardson number, J, for R — 3. The corresponding value of the Richardson number 
at the origin Ri(0) = RJ is indicated at the ordinate axis on the right. 




FIG. 15: (Color online) Comparison of the finite time Lyapunov exponent, In (<7 max (i))/t (solid), associated with the growth of 
optimal perturbations for optimizing times t = 50, t — 200 and t = 600, with the modal growth rate kd (dashed) for various 
Richardson numbers. The first three unstable Holmboe branches are shown. The HI branch is for Ri(0) < 1.42. The H2 
branch is for 10.46 < Ri(0) < 16.26; and H3 is for 31.04 < Ri(0) < 39.63. The optimal growth rate is almost constant for large 
Richardson numbers and is continuous across the islands of instability. Parameters are: k = 1 and R = 3. 
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0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

k 

FIG. 16: (Color online) Modal growth rate kd as a function of wavenumber k for center Richardson number Ri(0) = 5. The 
instability for 2.6 < k < 3.9 corresponds to the HI branch. The dots indicate the growth rate at k — 1, 2, 2.5, 4 and at k — 3.5 
when the maximum growth rate occurs. 



(a) k =1.75 (b) k =3.50 




FIG. 17: (Color online) Contours of the logarithm of the positive real part of the perturbation 5ft(i/>(z, i)e l x ) at z = 0.7 in 
the [x, t) plane for an initial perturbation in the form of the adjoint in the energy norm of the most unstable Holmboe mode 
at k = 3.5 for R = 3 and Ri(0) = 5. Panel (a): evolution under the dynamics with k = 1.75 when the flow is stable. Panel 
(b): evolution under the dynamics with k = 3.5 when the flow is unstable and this adjoint excites optimally the prograde 
Holmboe mode. The same propagation characteristics emerge at other wavenumbers for which the flow is stable. These 
Hovmoller diagrams indicate that when the flow is stable a quasi-mode emerges propagating with phase speed close to that of 
the Holmboe mode. 
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FIG. 18: (Color online) Energy evolution of a unit energy initial condition with the structure of the adjoint in the energy 
norm of the unstable HI mode at k — 3.5, R = 3 and Ri(0) = 5 with the dynamical operator (a) for the unstable k = 3.5 
(solid-squares), (b) for k = 1.0 (solid crosses), k = 1.75 (solid-dots), k = 2.0 (dashed), k — 2.5 (dash-dot) and k = 4.0 (solid) 
for which the flow is stable. When the flow is stable large amplitude propagating quasi-modes emerge as shown in FIG. 17. 
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FIG. 19: (Color online) Amplitude, |a(c r )|, of the absolute value of the expansion coefficients of the quasi-mode at k = 1.75 as 
a function of the phase speed of the modes of the perturbation operator at k = 1.75. All the modes of the perturbation operator 
are stable. The quasi-mode is excited by a unit energy initial condition with the structure of the adjoint of the unstable mode 
for k = 3.5. Notice the sharp peak at c r ~ 0.4278 which indicates the formation of long-lived quasi-modes that propagate with 
this phase speed (see FIG. 17(a)). For R = 3 and Ri(0) = 5. 
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FIG. 20: (Color online) The maximum energy response, ||R(u;)|| 2 , to harmonic forcing as a function of u/k for perturbations 
with k = 1, k = 2.5, k = 3.5 and k — 4.0. The dynamics have been rendered stable by addition of Rayleigh friction. The 
response is maximized at the phase speed of the quas-mode. Also shown is the resonant frequency response that would obtain 
if the operators were normal (dashed). The amplified response in energy is a measure of the non-normality of the operator. 
For stratification Ri(0) = 5 and R = 3. 
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FIG. 21: (Color online) Contours of modal growth rate kci as a function of wavenumber, k, and stratification parameter, J, 
of perturbations in Holmboe's idealized piecewise constant background state. At larger values of k and J there is a band of 
instability associated with the Holmboe instability. 
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FIG. 22: (Color online) Optimal energy growth that can be produced by the edge waves as a function of time for stratification 
parameter J = 1.5 for k = 4.2528, 4.24, 4. All k < 4.2528 are modally stable (fccj = 0) and k = 4.2528 is located at the 
stability boundary. The non-modal growth produced by the edge waves is substantial as the stability boundary is approached. 



